Systems and methods for modeling wellbore trajectories

ABSTRACT

Systems and methods for modeling wellbore trajectories, which can be used to model corresponding drillstring trajectories and transform the torque-drag drill string model into a full stiff-string formulation.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.12/337,408, which is hereby incorporated by reference, and claims thepriority of U.S. patent application Ser. No. 61/014,362, filed on Dec.17, 2007, which is incorporated herein by reference.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

Not applicable.

FIELD OF THE INVENTION

The present invention generally relates to modeling wellboretrajectories. More particularly, the present invention relates to theuse of spline functions, derived from drill string solutions, to modelwellbore trajectories.

BACKGROUND OF THE INVENTION

Wellbore trajectory models are used for two distinct purposes. The firstuse is planning the well location, which consists of determiningkick-off points, build and drop rates, and straight sections needed toreach a specified target. The second use is to integrate measuredinclination and azimuth angles to determine a well's location.

Various trajectory models have been proposed, with varying degrees ofsmoothness. The simplest model, the tangential model, consists ofstraight line sections. Thus, the slope of this model is discontinuousat survey points. The most commonly used model is the minimum curvaturemodel, which consists of circular-arc sections. This model hascontinuous slope, but discontinuous curvature. In fact, the minimumcurvature model argues that a wellbore would not necessarily havecontinuous curvature.

Analysis of drillstring loads is typically done with drillstringcomputer models. By far the most common method for drillstring analysisis the “torque-drag” model originally described in the Society ofPetroleum Engineers article “Torque and Drag in DirectionalWells—Prediction and Measurement” by Johancsik, C. A., Dawson, R. andFriesen, D. B., which was later translated into differential equationform as described in the article “Designing Well Paths to Reduce Dragand Torque” by Sheppard, M. C., Wick, C. and Burgess, T. M.

Torque-drag modeling refers to the torque and drag related todrillstring operation. Drag is the excess load compared to rotatingdrillstring weight, which may be either positive when pulling thedrillstring or negative while sliding into the well. This drag force isattributed to friction generated by drillstring contact with thewellbore. When rotating, this same friction will reduce the surfacetorque transmitted to the bit. Being able to estimate the frictionforces is useful when planning a well or analysis afterwards. Because ofthe simplicity and general availability of the torque-drag model, it hasbeen used extensively for planning and in the field. Field experienceindicates that this model generally gives good results for many wells,but sometimes performs poorly.

In the standard torque-drag model, the drillstring trajectory is assumedto be the same as the wellbore trajectory, which is a reasonableassumption considering that surveys are taken within the drillstring.Contact with the wellbore is assumed to be continuous. However, giventhat the most common method for determining the wellbore trajectory isthe minimum curvature method, the wellbore shape is less than idealbecause the bending moment is not continuous and smooth at surveypoints. This problem is dealt with by neglecting bending moment but, asa result of this assumption, some of the contact force is alsoneglected.

Therefore, there is a need for a new wellbore trajectory model that hassufficient smoothness to model the drillstring trajectory.

There is a further need to provide a new wellbore trajectory model thattransforms the simple torque-drag drill string model into a fullstiff-string formulation because, in this formulation, drill stringbending and shear forces arise that cannot be determined correctly withconventional wellbore trajectory models.

SUMMARY OF THE INVENTION

The present invention meets the above needs and overcomes one or moredeficiencies in the prior art by providing systems and methods formodeling a wellbore trajectory, which can be used to model thecorresponding drillstring trajectory and transform the torque-drag drillstring model into a full stiff-string formulation.

In one embodiment, the present invention includes a computer implementedmethod for modeling a wellbore trajectory, which comprises: i)calculating a tangent vector interpolation function for each intervalbetween two or more survey points within a wellbore; and (ii)determining the wellbore trajectory using a computer processor and eachtangent vector interpolation function in a torque-drag drillstringmodel.

In another embodiment, the present invention includes a non-transitorycomputer readable medium having computer executable instructions formodeling a wellbore trajectory. The instructions are executable toimplement: i) calculating a tangent vector interpolation function foreach interval between two or more survey points within a wellbore; and(ii) determining the wellbore trajectory using each tangent vectorinterpolation function in a torque-drag drillstring model.

Additional aspects, advantages and embodiments of the invention willbecome apparent to those skilled in the art from the followingdescription of the various embodiments and related drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention is described below with references to theaccompanying drawings in which like elements are referenced with likereference numerals, and in which:

FIG. 1 is a block diagram illustrating one embodiment of a system forimplementing the present invention.

FIG. 2 is a graphical illustration comparing the analytic model, theminimum curvature model and the spline model of the present inventionfor a circular-arc wellbore trajectory.

FIG. 3 is a graphical illustration comparing the analytic model, theminimum curvature model and the spline model of the present inventionfor a catenary wellbore trajectory.

FIG. 4 is a graphical illustration comparing the analytic model, theminimum curvature model and the spline model of the present inventionfor a helix wellbore trajectory.

FIG. 5 is a graphical illustration comparing the rate-of-change ofcurvature between an analytic model and the spline model of the presentinvention for a catenary wellbore trajectory.

FIG. 6 is a graphical illustration comparing the torsion between ananalytic model and the spline model of the present invention for a helixwellbore trajectory.

FIG. 7 illustrates the test case wellbore used in Example 1.

FIG. 8 is a graphical illustration comparing the bending moment betweenthe minimum curvature model and the spline model of the presentinvention for the test case wellbore used in Example 1.

FIG. 9A is a graphical illustration (vertical view) of the short radiuswellpath used in Example 2.

FIG. 9B is a graphical illustration (North/East view) of the shortradius wellpath used in Example 2.

FIG. 10 is a graphical illustration comparing the short radius contactforce between a constant curvature model and the spline model of thepresent invention for the wellpath used in Example 2.

FIG. 11 is a graphical illustration comparing the short radius bendingmoment between a constant curvature model and the spline model of thepresent invention for the wellpath used in Example 2.

FIG. 12 is a flow diagram illustrating one embodiment of a method forimplementing the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The subject matter of the present invention is described withspecificity, however, the description itself is not intended to limitthe scope of the invention. The subject matter thus, might also beembodied in other ways, to include different steps or combinations ofsteps similar to the ones described herein, in conjunction with otherpresent or future technologies. Moreover, although the term “step” maybe used herein to describe different elements of methods employed, theterm should not be interpreted as implying any particular order among orbetween various steps herein disclosed unless otherwise expresslylimited by the description to a particular order.

System Description

The present invention may be implemented through a computer-executableprogram of instructions, such as program modules, generally referred toas software applications or application programs executed by a computer.The software may include, for example, routines, programs, objects,components, and data structures that perform particular tasks orimplement particular abstract data types. The software forms aninterface to allow a computer to react according to a source of input.WELLPLAN™, which is a commercial software application marketed byLandmark Graphics Corporation, may be used as an interface applicationto implement the present invention. The software may also cooperate withother code segments to initiate a variety of tasks in response to datareceived in conjunction with the source of the received data. Thesoftware may be stored and/or carried on any variety of memory mediasuch as CD-ROM, magnetic disk, bubble memory and semiconductor memory(e.g., various types of RAM or ROM). Furthermore, the software and itsresults may be transmitted over a variety of carrier media such asoptical fiber, metallic wire, free space and/or through any of a varietyof networks such as the Internet.

Moreover, those skilled in the art will appreciate that the inventionmay be practiced with a variety of computer-system configurations,including hand-held devices, multiprocessor systems,microprocessor-based or programmable-consumer electronics,minicomputers, mainframe computers, and the like. Any number ofcomputer-systems and computer networks are acceptable for use with thepresent invention. The invention may be practiced indistributed-computing environments where tasks are performed byremote-processing devices that are linked through a communicationsnetwork. In a distributed-computing environment, program modules may belocated in both local and remote computer-storage media including memorystorage devices. The present invention may therefore, be implemented inconnection with various hardware, software or a combination thereof, ina computer system or other processing system.

Referring now to FIG. 1, a block diagram of a system for implementingthe present invention on a computer is illustrated. The system includesa computing unit, sometimes referred to as a computing system, whichcontains memory, application programs, a client interface, and aprocessing unit. The computing unit is only one example of a suitablecomputing environment and is not intended to suggest any limitation asto the scope of use or functionality of the invention.

The memory primarily stores the application programs, which may also bedescribed as program modules containing computer-executableinstructions, executed by the computing unit for implementing themethods described herein and illustrated in FIGS. 2-12. The memorytherefore, includes a Wellbore Trajectory Module, which enables themethods illustrated and described in reference to FIGS. 2-12, andWELLPLAN™.

Although the computing unit is shown as having a generalized memory, thecomputing unit typically includes a variety of computer readable media.By way of example, and not limitation, computer readable media maycomprise computer storage media and communication media. The computingsystem memory may include computer storage media in the form of volatileand/or nonvolatile memory such as a read only memory (ROM) and randomaccess memory (RAM). A basic input/output system (BIOS), containing thebasic routines that help to transfer information between elements withinthe computing unit, such as during start-up, is typically stored in ROM.The RAM typically contains data and/or program modules that areimmediately accessible to, and/or presently being operated on by, theprocessing unit. By way of example, and not limitation, the computingunit includes an operating system, application programs, other programmodules, and program data.

The components shown in the memory may also be included in otherremovable/nonremovable, volatile/nonvolatile computer storage media. Forexample only, a hard disk drive may read from or write to nonremovable,nonvolatile magnetic media, a magnetic disk drive may read from or writeto a removable, non-volatile magnetic disk, and an optical disk drivemay read from or write to a removable, nonvolatile optical disk such asa CD ROM or other optical media. Other removable/non-removable,volatile/non-volatile computer storage media that can be used in theexemplary operating environment may include, but are not limited to,magnetic tape cassettes, flash memory cards, digital versatile disks,digital video tape, solid state RAM, solid state ROM, and the like. Thedrives and their associated computer storage media discussed abovetherefore, store and/or carry computer readable instructions, datastructures, program modules and other data for the computing unit.

A client may enter commands and information into the computing unitthrough the client interface, which may be input devices such as akeyboard and pointing device, commonly referred to as a mouse, trackballor touch pad. Input devices may include a microphone, joystick,satellite dish, scanner, or the like.

These and other input devices are often connected to the processing unitthrough the client interface that is coupled to a system bus, but may beconnected by other interface and bus structures, such as a parallel portor a universal serial bus (USB). A monitor or other type of displaydevice may be connected to the system bus via an interface, such as avideo interface. In addition to the monitor, computers may also includeother peripheral output devices such as speakers and printer, which maybe connected through an output peripheral interface.

Although many other internal components of the computing unit are notshown, those of ordinary skill in the art will appreciate that suchcomponents and their interconnection are well known.

Method Description

Unlike prior wellbore trajectory models, the present invention proceedsfrom the concept that the trajectory given by the survey measurementsmade within the drillstring is the trajectory of the drillstring, whichmust have continuity of bending moment proportional to curvature. Thenomenclature used herein is described in the Society of PetroleumEngineers article “Drillstring Solutions Improve the Torque-Drag Model”by Mitchell, Robert F. (“SPE 112623”), which is incorporated herein byreference and repeated in Table 1 below.

TABLE 1 {right arrow over (b)} binormal vector {tilde over (b)} specialbinormal vector E Young's elastic modulus (psf) F the effective axialforce (lbf.) {tilde over (F)} F − EIκ² I moment of inertia (ft⁴)

unit vector in east direction

unit vector in north direction

unit vector in downward direction {right arrow over (n)} normal vector ñspecial normal vector s measured depth (ft) {right arrow over (t)}tangent vector

spline tangent vector function

position vector, (ft)

initial position vector, increment j (ft) α_(j) coefficient in ñdirection (ft⁻¹) β_(j) coefficient in {tilde over (b)} direction (ft⁻¹)Δs_(j) s_(j+1) − s_(j) (ft) λ_(j) Coefficient in spline functions ε_(j)angle between {right arrow over (n)} and ñ κ wellbore curvature (ft⁻¹) φwellbore trajectory inclination angle θ wellbore trajectory azimuthangle ξ_(j) (s − s_(j))/(s_(j+1) − s_(j)) ′ d/ds ^(iv) d⁴/ds⁴ subscriptsj survey point

The use of cubic splines is well known in the art for achieving highercontinuity in a trajectory model. If, for example, a table of {x_(i),y_(i)} is used, intermediate values of y as a function of x may bedetermined by linear interpolation:

$\begin{matrix}{{y(x)} = {{y_{j}\left( \frac{x_{j + 1} - x}{x_{j + 1} - x_{j}} \right)} + {y_{j + 1}\left( \frac{x - x_{j}}{x_{j + 1} - x_{j}} \right)}}} & (1)\end{matrix}$where the interpolation occurs between x_(j) and x_(j+1). If it isdesired that the interpolation have smooth first and second derivativesat the x_(j) points, the interpolation may be:y(x)=y _(j) f ₁(x)+y″ _(j) f ₂(x)+y _(j+1) f ₃(x)+y″ _(j+1) f ₄(x)  (2)where the functions f_(j) are devised so that:y(x _(j))=y _(j)y(x _(j+1))=y _(j+1)y″(x _(j))=y″ _(j)y″(x _(j+1))=y″ _(j+1)  (3)

In the classic cubic spline formulation, the f_(j) are cubic functionsof x and the unknown coefficients y″_(j) are determined by requiringcontinuity of the first derivatives of y(x) at each x_(j). Here thefunctions in equation (2) need not be cubic functions. They must onlysatisfy equations (3). The use of spline formulations such as, forexample, cubic splines and tangent splines to model wellboretrajectories is well known in the art. The determination of the wellboretrajectory from survey data, however, is not. Furthermore, the use ofconventional splines, as applied to a three-dimensional curve, will notsatisfy equation (5) and equation (6).

Once survey data is obtained, the tangent vector {right arrow over(t)}_(j) at each survey point j can be calculated. One formula forinterpolating the tangent vectors is:

$\begin{matrix}{{{\overset{\rightharpoonup}{t}}_{j}(s)} = {\quad{{\begin{matrix}{{\overset{->}{T}}_{j}(s)} \\\sqrt{{{\overset{\rightharpoonup}{T}}_{j}(s)} \cdot {{\overset{\rightharpoonup}{T}}_{j}(s)}}\end{matrix}{T_{j}(s)}} = {{{\overset{\rightharpoonup}{t}}_{j}{f_{1j}(s)}} + {\kappa_{j}{\overset{\rightharpoonup}{n}}_{j}{f_{2j}(s)}} + {{\overset{\rightharpoonup}{t}}_{j + 1}{f_{3j}(s)}} + {\kappa_{j + 1}{\overset{\rightharpoonup}{n}}_{j + 1}{f_{4j}(s)}}}}}} & (4)\end{matrix}$where s is measured depth, κ_(j) is the curvature at s_(j), and {rightarrow over (n)}_(j) is the normal vector at s_(j). This formulation hastwo purposes. The first purpose is to satisfy the Frenet equation for acurve (by suitable choice of functions f_(ij)):

$\begin{matrix}{\frac{\mathbb{d}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s} = {{\kappa(s)}{\overset{\rightharpoonup}{n}(s)}}} & (5)\end{matrix}$The second reason is to insure that s is indeed measured depth. Thisrequirement means:du ₁ ² +du ₂ ² +du ₃ ² =ds ²(an incremental change of position equals the incremental arc length)or, in terms of the tangent vectors:

$\begin{matrix}{{\left( \frac{\mathbb{d}u_{1}}{\mathbb{d}s} \right)^{2} + \left( \frac{\mathbb{d}u_{2}}{\mathbb{d}s} \right)^{2} + \left( \frac{\mathbb{d}u_{3}}{\mathbb{d}s} \right)^{2}} = {{\frac{\mathbb{d}\overset{\rightharpoonup}{u}}{\mathbb{d}s} \cdot \frac{\mathbb{d}\overset{\rightharpoonup}{u}}{\mathbb{d}s}} = {{\overset{\rightharpoonup}{t} \cdot \overset{\rightharpoonup}{t}} = 1}}} & (6)\end{matrix}$

As demonstrated in the following section, equation (4) satisfies thiscondition. The details for determining the unknowns in equation (4),which are the normal vectors and the curvatures, are also addressed inthe following section.

Spline Wellbore Trajectory

The normal method for determining the well path is to use some type ofsurveying instrument to measure the inclination and azimuth at variousdepths and then to calculate the trajectory. At each survey point j,inclination angle φ_(j) and azimuth angle θ_(j) are measured, as well asthe course length Δs=s_(j)=s_(j+1)−s_(j) between survey points. Eachsurvey point j therefore, includes survey data comprising an inclinationangle φ_(j), an azimuth angle θ_(j) and a measured depth s. These angleshave been corrected (i) to true north for a magnetic survey or (ii) fordrift if a gyroscopic survey. The survey angles define the tangent{right arrow over (t)}_(j) to the trajectory at each survey point jwhere the tangent vector is defined in terms of inclination φ_(j) andazimuth θ_(j) in the following formulas:{right arrow over (t)} _(j) ●{right arrow over (i)}_(N)=cos(θ_(j))sin(φ_(j)){right arrow over (t)} _(j) ●{right arrow over (i)}_(E)=sin(θ_(j))sin(φ_(j)){right arrow over (t)} _(j) ●{right arrow over (i)} _(z)=cos(φ_(j))  (7)

If it was known how the angles φ and θ varied between survey points, orequivalently, if it was known how the tangent vectors varied betweensurvey points, then the trajectory could be determined by integratingthe tangent vector:

$\begin{matrix}{{{{\overset{\rightharpoonup}{t}}_{j} = \frac{\mathbb{d}{\overset{\rightharpoonup}{u}}_{j}}{\mathbb{d}s}},{so}}{{{\overset{\rightharpoonup}{u}}_{j}(s)} = {{\overset{\rightharpoonup}{u}\frac{o}{j}} + {\int_{j}{{\overset{\rightharpoonup}{t}}_{j}{\mathbb{d}s}}}}}} & (8)\end{matrix}$

Given tangent vectors and {right arrow over (t)}_(j) and {right arrowover (t)}_(j+1) and associated normal vectors {right arrow over (n)}_(j)and {right arrow over (n)}_(j+1), a tangent vector interpolationfunction connecting these vectors can be created. First, a set ofinterpolation functions f_(ij)(s), s in [s_(j), s_(j+1)], with thefollowing properties, will be needed:

$\begin{matrix}{{{{f_{1j}\left( s_{j} \right)} = 1},{\frac{\mathbb{d}{f_{1j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{1j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{2j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{2j}\left( s_{j} \right)}}{\mathbb{d}s} = 1},{{f_{2j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{3j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{3j}\left( s_{j + 1} \right)} = 1},{\frac{\mathbb{d}{f_{3j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{4j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{4j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{4j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 1}}} & (9)\end{matrix}$

There are a variety of functions that satisfy equations (9). If thespline function T_(j)(ξ) is defined as:

$\begin{matrix}{{T_{j}(\xi)} = {{{\overset{\rightharpoonup}{t}}_{j}{f_{1j}(s)}} + {\kappa_{j}{\overset{\rightharpoonup}{n}}_{j}{f_{2j}(s)}} + {{\overset{\rightharpoonup}{t}}_{j + 1}{f_{3j}(s)}} + {\kappa_{j + 1}{\overset{\rightharpoonup}{n}}_{j + 1}{f_{4j}(s)}}}} & (10)\end{matrix}$it becomes clear that:

$\begin{matrix}{{{{\overset{\rightharpoonup}{T}}_{j}\left( s_{j} \right)} = {\overset{->}{t}}_{j}}{{{\overset{\rightharpoonup}{T}}_{j}\left( s_{j + 1} \right)} = {\overset{\rightharpoonup}{t}}_{j + 1}}{{\frac{\mathbb{d}{\overset{\rightharpoonup}{T}}_{j}}{\mathbb{d}s}\left( s_{j} \right)} = {\kappa_{j}{\overset{\rightharpoonup}{n}}_{j}}}{{\frac{\mathbb{d}{\overset{\rightharpoonup}{T}}_{j}}{\mathbb{d}s}\left( s_{j + 1} \right)} = {\kappa_{j + 1}{\overset{\rightharpoonup}{n}}_{j + 1}}}} & (11)\end{matrix}$The function T_(j) satisfies the Frenet equation:

$\begin{matrix}{\frac{\mathbb{d}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s} = {{\kappa(s)}{\overset{\rightharpoonup}{n}(s)}}} & (12)\end{matrix}$for a tangent vector at s=s_(j) and s_(j+1). However, T_(j) is not atangent vector because it is not a unit vector. This can be corrected bynormalizing Tj:

$\begin{matrix}{{{\overset{\rightharpoonup}{t}}_{j}(s)} = {\quad\begin{matrix}{{\overset{->}{T}}_{j}(s)} \\\sqrt{{{\overset{\rightharpoonup}{T}}_{j}(s)} \cdot {{\overset{\rightharpoonup}{T}}_{j}(s)}}\end{matrix}}} & (13)\end{matrix}$where it is shown that equation (12) is still satisfied. In order toevaluate the curvatures κ_(j), equation (13) is differentiated twice andevaluated at s=s_(j) and s_(j+1):

$\begin{matrix}{\mspace{79mu}{{{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{t}}_{j}} = {- \kappa_{j}^{2}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{n}}_{j}} = {{\kappa_{j}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\rightharpoonup}{n}}_{j} \cdot {\overset{\rightharpoonup}{t}}_{j + 1}}\frac{\mathbb{d}^{2}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {\kappa_{j + 1}{n_{j + 1} \cdot {\overset{\rightharpoonup}{n}}_{j}}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{b}}_{j}} = {{{{\overset{\rightharpoonup}{t}}_{j + 1} \cdot {\overset{\rightharpoonup}{b}}_{j}}\frac{\mathbb{d}^{2}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {\kappa_{j + 1}{{\overset{\rightharpoonup}{n}}_{j + 1} \cdot {\overset{\rightharpoonup}{b}}_{j}}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}}}}}\mspace{79mu}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{t}}_{j + 1}} = {- \kappa_{j + 1}^{2}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{n}}_{j + 1}} = {{\kappa_{j + 1}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\rightharpoonup}{n}}_{j + 1} \cdot {\overset{\rightharpoonup}{t}}_{j}}\frac{\mathbb{d}^{2}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {\kappa_{j}{{\overset{\rightharpoonup}{n}}_{j} \cdot n_{j + 1}}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{b}}_{j + 1}} = {{{{\overset{\rightharpoonup}{t}}_{j} \cdot {\overset{\rightharpoonup}{b}}_{j + 1}}\frac{\mathbb{d}^{2}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {\kappa_{j}{{\overset{\rightharpoonup}{n}}_{j} \cdot {\overset{\rightharpoonup}{b}}_{j + 1}}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}}}}}} & (14)\end{matrix}$Using the Frenet equation (12) and

$\begin{matrix}{{\frac{\mathbb{d}{\overset{\rightharpoonup}{n}(s)}}{\mathbb{d}s} = {{{- {\kappa(s)}}{\overset{\rightharpoonup}{t}(s)}} + {{\tau(s)}{\overset{\rightharpoonup}{b}(s)}}}}{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s^{2}} = {{{- {\kappa^{2}(s)}}{\overset{\rightharpoonup}{t}(s)}} + {{\kappa^{\prime}(s)}{\overset{\rightharpoonup}{n}(s)}} + {{\kappa(s)}{\tau(s)}{\overset{\rightharpoonup}{b}(s)}}}}} & (15)\end{matrix}$it is evident that:

$\begin{matrix}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{t}}_{j}} = {- \kappa_{j}^{2}}} & \left( {16a} \right) \\{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{n}}_{j}} = \frac{\mathbb{d}\kappa_{j}}{\mathbb{d}s}} & \left( {16b} \right) \\{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{b}}_{j}} = {\kappa_{j}\tau_{j}}} & \left( {16c} \right) \\{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{t}}_{j + 1}} = {- \kappa_{j + 1}^{2}}} & \left( {16d} \right) \\{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{n}}_{j + 1}} = \frac{\mathbb{d}\kappa_{j + 1}}{\mathbb{d}s}} & \left( {16e} \right) \\{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\rightharpoonup}{b}}_{j + 1}} = {\kappa_{j + 1}\tau_{j + 1}}} & \left( {16f} \right)\end{matrix}$The Frenet formulae, equation (15), are identically satisfied byequation (16a) and equation (16d). Before this set of equations can besolved for curvatures κ_(j), a representation for the normal vector({right arrow over (n)}_(j)) and the binormal vector ({right arrow over(b)}_(j)) is needed. The tangent vector is defined by the inclinationangle (φ_(j)) and the azimuth angle (θ_(j)) in the following way:

$\begin{matrix}{{\overset{->}{t}}_{j} = \begin{bmatrix}{\sin\;\varphi_{j}\cos\;\vartheta_{j}} \\{\sin\;\varphi_{j}\sin\;\vartheta_{j}} \\{\cos\;\varphi_{j}}\end{bmatrix}} & (17)\end{matrix}$Then the Frenet equation (7) requires:

$\begin{matrix}\begin{matrix}{{\frac{\mathbb{d}}{\mathbb{d}s}{\overset{\rightharpoonup}{t}}_{j}} = {{\begin{bmatrix}{\cos\;\varphi_{j}\cos\;\vartheta_{j}} \\{\cos\;\varphi_{j}\sin\;\vartheta_{j}} \\{{- \sin}\;\varphi_{j}}\end{bmatrix}\frac{\mathbb{d}}{\mathbb{d}s}\varphi_{j}} + {\begin{bmatrix}{{- \sin}\;\vartheta_{j}} \\{\cos\;\vartheta_{j}} \\0\end{bmatrix}\sin\;\varphi_{j}\frac{\mathbb{d}}{\mathbb{d}s}\vartheta_{j}}}} \\{= {\kappa_{j}{\overset{\rightharpoonup}{n}}_{j}}}\end{matrix} & (18)\end{matrix}$From equation (12), the equation for the curvature κ_(j) becomes:

$\begin{matrix}{\kappa_{j} = \sqrt{\left( {\frac{\mathbb{d}}{\mathbb{d}s}\varphi_{j}} \right)^{2} + {\sin^{2}\;{\varphi_{j}\left( {\frac{\mathbb{d}}{\mathbb{d}s}\vartheta_{j}} \right)}^{2}}}} & (19)\end{matrix}$We define the following quantities found in equation (18):

$\begin{matrix}{{{\overset{\sim}{n}}_{j} = \begin{bmatrix}{\cos\;\varphi_{j}\cos\;\vartheta_{j}} \\{\cos\;\varphi_{j}\sin\;\vartheta_{j}} \\{{- \sin}\;\varphi_{j}}\end{bmatrix}}{{\overset{\sim}{b}}_{j} = \begin{bmatrix}{{- \sin}\;\vartheta_{j}} \\{\cos\;\vartheta_{j}} \\0\end{bmatrix}}} & (20)\end{matrix}$These vectors are useful in defining the normal and binormal vectors.

As provided above, {right arrow over (t)}_(j), ñ_(j), and {tilde over(b)}_(j) form a right-handed coordinate system at s_(j). The normalvector ({right arrow over (n)}_(j)) and the binormal vector ({rightarrow over (b)}_(j)) can be defined by rotation through the angle ε_(j)around the tangent vector:{right arrow over (n)} _(j) =ñ _(j) cos ε_(j) +{tilde over (b)} _(j) sinε_(j){right arrow over (b)} _(j) =−ñ _(j) sin ε_(j) +{tilde over (b)} _(j)cos ε_(j)  (21)Then {right arrow over (n)}_(j) is a unit vector consistent with Frenetequation (5), given:

$\begin{matrix}{{\cos\; ɛ_{j}} = {{\frac{1}{\kappa_{j}}\frac{\mathbb{d}}{\mathbb{d}s}\varphi_{j}\mspace{14mu}{and}\mspace{14mu}\sin\mspace{14mu} ɛ_{j}} = {\frac{\sin\;\varphi_{j}}{\kappa_{j}}\frac{\mathbb{d}}{\mathbb{d}s}\vartheta_{j}}}} & (22)\end{matrix}$The variables κ_(j) and ε_(j) are not the most convenient choicesbecause of the nonlinearity introduced by the sine and cosine functions.An alternate selection may be:

$\begin{matrix}{{{\kappa_{j}{\overset{\rightharpoonup}{n}}_{j}} = {{\alpha_{j}{\overset{\sim}{n}}_{j}} + {\beta_{j}{\overset{\sim}{b}}_{j}}}}{\alpha_{j} = {\kappa_{j}\cos\; ɛ_{j}}}{\beta_{j} = {\kappa_{j}\sin\; ɛ_{j}}}{\kappa_{j} = \sqrt{\alpha_{j}^{2} + \beta_{j}^{2}}}{ɛ_{j} = {\tan^{- 1}\frac{\beta_{j}}{\alpha_{j}}}}} & (23)\end{matrix}$Equations (16a)-(16f) can be rewritten in terms of the vectors ñ and{tilde over (b)} to give:

$\begin{matrix}{{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\sim}{n}}_{j}} = {{\alpha_{j}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{n}}_{j} \cdot {\overset{\rightharpoonup}{t}}_{j + 1}}\frac{\mathbb{d}^{2}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{n}}_{j} \cdot \left\lbrack {{\alpha_{j + 1}{\overset{\sim}{n}}_{j + 1}} + {\beta_{j + 1}{\overset{\sim}{b}}_{j + 1}}} \right\rbrack}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\sim}{b}}_{j}} = {{\beta_{j}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{b}}_{j} \cdot {\overset{\rightharpoonup}{t}}_{j + 1}}\frac{\mathbb{d}^{2}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{b}}_{j} \cdot \left\lbrack {{\alpha_{j + 1}{\overset{\sim}{n}}_{j + 1}} + {\beta_{j + 1}{\overset{\sim}{b}}_{j + 1}}} \right\rbrack}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s^{2}}}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\sim}{n}}_{j + 1}} = {{\alpha_{j + 1}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{n}}_{j + 1} \cdot {\overset{\rightharpoonup}{t}}_{j}}\frac{\mathbb{d}^{2}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{n}}_{j + 1} \cdot \left( {{\alpha_{j}{\overset{\sim}{n}}_{j}} + {\beta_{j}{\overset{\sim}{b}}_{j}}} \right)}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}}}}{{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}} \cdot {\overset{\sim}{b}}_{j + 1}} = {{\beta_{j + 1}\frac{\mathbb{d}^{2}{f_{4j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{b}}_{j + 1} \cdot {\overset{\rightharpoonup}{t}}_{j}}\frac{\mathbb{d}^{2}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}} + {{{\overset{\sim}{b}}_{j + 1} \cdot \left( {{\alpha_{j}{\overset{\sim}{n}}_{j}} + {\beta_{j}{\overset{\sim}{b}}_{j}}} \right)}\frac{\mathbb{d}^{2}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s^{2}}}}}} & (24)\end{matrix}$Continuity of d²{right arrow over (t)}/ds² at survey points requires forj=2, N−1:

$\begin{matrix}{{{\left\lbrack {{\alpha_{j - 1}{{\overset{\sim}{n}}_{j} \cdot {\overset{\sim}{n}}_{j - 1}}} + {\beta_{j - 1}{{\overset{\sim}{n}}_{j} \cdot {\overset{\sim}{b}}_{j - 1}}}} \right\rbrack\frac{\mathbb{d}^{2}f_{{2j} - 1}}{\mathbb{d}s^{2}}} + {\alpha_{j}\left( {\frac{\mathbb{d}^{2}f_{2j}}{\mathbb{d}s^{2}} + \frac{\mathbb{d}^{2}f_{{4j} - 1}}{\mathbb{d}s^{2}}} \right)} + {\left\lbrack {{\alpha_{j + 1}{{\overset{\sim}{n}}_{j} \cdot {\overset{\sim}{n}}_{j + 1}}} + {\beta_{j + 1}{{\overset{\sim}{n}}_{j} \cdot {\overset{\sim}{b}}_{j + 1}}}} \right\rbrack\frac{\mathbb{d}^{2}f_{4j}}{\mathbb{d}s^{2}}}} = {{{{{\overset{\sim}{n}}_{j} \cdot {\left( {{{\overset{\rightharpoonup}{t}}_{j + 1}\frac{\mathbb{d}^{2}f_{3j}}{\mathbb{d}s^{2}}} - {{\overset{\rightharpoonup}{t}}_{j - 1}\frac{\mathbb{d}^{2}f_{{1j} - 1}}{\mathbb{d}s^{2}}}} \right)\left\lbrack {{\alpha_{j - 1}{{\overset{\sim}{b}}_{j} \cdot {\overset{\sim}{n}}_{j - 1}}} + {\beta_{j - 1}{{\overset{\sim}{b}}_{j} \cdot {\overset{\sim}{b}}_{j - 1}}}} \right\rbrack}}\frac{\mathbb{d}^{2}f_{{2j} - 1}}{\mathbb{d}s^{2}}} + {\beta_{j}\left( {\frac{\mathbb{d}^{2}f_{2j}}{\mathbb{d}s^{2}} + \frac{\mathbb{d}^{2}f_{{4j} - 1}}{\mathbb{d}s^{2}}} \right)} + {\left\lbrack {{\alpha_{j + 1}{{\overset{\sim}{b}}_{j} \cdot {\overset{\sim}{n}}_{j + 1}}} + {\beta_{j + 1}{{\overset{\sim}{b}}_{j} \cdot {\overset{\sim}{b}}_{j + 1}}}} \right\rbrack\frac{\mathbb{d}^{2}f_{4j}}{\mathbb{d}s^{2}}}} = {{\overset{\sim}{b}}_{j} \cdot \left( {{{\overset{\rightharpoonup}{t}}_{j + 1}\frac{\mathbb{d}^{2}f_{3j}}{\mathbb{d}s^{2}}} - {{\overset{\rightharpoonup}{t}}_{j - 1}\frac{\mathbb{d}^{2}f_{{1j} - 1}}{\mathbb{d}s^{2}}}} \right)}}} & (25)\end{matrix}$The set of equations (25) together with boundary conditions defined atthe initial and end points form a diagonally dominant block tridiagonalset of equations that are relatively easy to solve. Notably, by alsosolving for α_(j) and β_(j), the system has also solved for dφ_(j)/dsand dθ_(j)/ds through equation (23). Further, there is no ambiguityabout the magnitude of θ_(j) (±nπ) in the definition of thesederivatives.

There is therefore, a need for expressions for the parameters κ, τ, andκ′ that appear in the torque-drag equilibrium equations.

Recalling the Frenet formulae (equations (12) and (15)):

$\begin{matrix}{{\frac{\mathbb{d}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s} = {{\kappa(s)}{\overset{\rightharpoonup}{n}(s)}}}{\frac{\mathbb{d}{\overset{\rightharpoonup}{n}(s)}}{\mathbb{d}s} = {{{- {\kappa(s)}}{\overset{\rightharpoonup}{t}(s)}} + {{\tau(s)}{\overset{\rightharpoonup}{b}(s)}}}}{\frac{\mathbb{d}^{2}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s^{2}} = {{{- {\kappa^{2}(s)}}{\overset{\rightharpoonup}{t}(s)}} + {{\kappa^{\prime}(s)}{\overset{\rightharpoonup}{n}(s)}} + {{\kappa(s)}{\tau(s)}{\overset{\rightharpoonup}{b}(s)}}}}{{{\overset{\rightharpoonup}{t}(s)} \times \frac{\mathbb{d}{\overset{\rightharpoonup}{t}(s)}}{\mathbb{d}s}} = {{{\overset{\rightharpoonup}{t}(s)} \times {\kappa(s)}{\overset{\rightharpoonup}{n}(s)}} = {{\kappa(s)}{\overset{\rightharpoonup}{b}(s)}}}}} & (26)\end{matrix}$it is determined that:

$\begin{matrix}{{{\kappa(s)} = \sqrt{\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{->}{t}}_{j}(s)} \cdot \frac{\mathbb{d}}{\mathbb{d}s}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}}{{{\kappa(s)}\frac{\mathbb{d}}{\mathbb{d}s}{\kappa(s)}} = {\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \frac{\mathbb{d}^{2}}{\mathbb{d}s^{2}}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}}{{{\kappa(s)}^{2}{\tau(s)}} = {\frac{\mathbb{d}^{2}}{\mathbb{d}s^{2}}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \left\lbrack {{{\overset{\rightharpoonup}{t}}_{j}(s)} \times \frac{\mathbb{d}}{\mathbb{d}s}{{\overset{\rightharpoonup}{t}}_{j}(\xi)}} \right\rbrack}}}} & (27)\end{matrix}$If κ is non-zero at a given point, then:

$\begin{matrix}{{{\kappa(s)} = \sqrt{\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \frac{\mathbb{d}}{\mathbb{d}s}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}}{{\frac{\mathbb{d}}{\mathbb{d}s}{\kappa(s)}} = \frac{\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \frac{\mathbb{d}^{2}}{\mathbb{d}s^{2}}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}{\sqrt{\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \frac{\mathbb{d}}{\mathbb{d}s}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}}}{{\tau(s)} = \frac{\frac{\mathbb{d}^{2}}{\mathbb{d}s^{2}}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \left\lbrack {{{\overset{\rightharpoonup}{t}}_{j}(s)} \times \frac{\mathbb{d}}{\mathbb{d}s}{{\overset{\rightharpoonup}{t}}_{j}(s)}} \right\rbrack}}{\frac{\mathbb{d}}{\mathbb{d}s}{{{\overset{\rightharpoonup}{t}}_{j}(s)} \cdot \frac{\mathbb{d}}{\mathbb{d}s}}{{\overset{\rightharpoonup}{t}}_{j}(s)}}}} & (28)\end{matrix}$

Since the system is intended to model drillstrings, the best choice forinterpolating functions (f_(ij)) are solutions to actual drillstringproblems.

The equation for the mechanical equilibrium of a weightless elastic rodwith large displacement is:EI{right arrow over (u)} ^(iv)−[(F−EIκ ²){right arrow over(u)}′]′={right arrow over (0)}  (29)where EI is the bending stiffness, F is the axial force (tensionpositive), and κ is the curvature of the rod. Looking at a smallinterval of s, F and κ are roughly constant, so the solution to equation(7) becomes:u(s)=c ₀ +c ₁ s+c ₂ sin h(λs)+c ₃ cos h(λs)when: EIλ ² =F−EIκ ²>0  (30a)u(s)=c ₀ +c ₁ s+c ₂ sin(λs)+c ₃ cos(λs)when: EIλ ² =EIκ ² −F>0  (30b)u(s)=c ₀ +c ₁ s+c ₂ s ² c ₃ s ³when: EIκ ² −F=0  (30c)where the c₀-c₃ are four constants to be determined. The third equationis a cubic equation, so cubic splines are a candidate solution, eventhough they represent a special case of zero axial loads. Equation (30a)can be used to define what are known as tension-splines and equation(30b) may be used to define “compression” splines. This is demonstratedin the following section using drillstring solutions as interpolationfunctions.

Drillstring Solutions as Interpolation Functions

As demonstrated in the Spline Wellbore Trajectory section above, a setof interpolation functions f_(ij)(s), s in [s_(j), s_(j+1)], is neededwith the following properties:

$\begin{matrix}{{{{f_{1j}\left( s_{j} \right)} = 1},{\frac{\mathbb{d}{f_{1j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{1j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{1j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{2j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{2j}\left( s_{j} \right)}}{\mathbb{d}s} = 1},{{f_{2j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{2j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{3j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{3j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{3j}\left( s_{j + 1} \right)} = 1},{\frac{\mathbb{d}{f_{3j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 0}}{{{f_{4j}\left( s_{j} \right)} = 0},{\frac{\mathbb{d}{f_{4j}\left( s_{j} \right)}}{\mathbb{d}s} = 0},{{f_{4j}\left( s_{j + 1} \right)} = 0},{\frac{\mathbb{d}{f_{4j}\left( s_{j + 1} \right)}}{\mathbb{d}s} = 1}}} & (31)\end{matrix}$For example, the following cubic functions satisfy the requirements ofequation (31):

$\begin{matrix}{{{f_{1j}(s)} = {1 + {\left( {{2\xi} - 3} \right)\xi^{2}}}}{{f_{2j}(s)} = {{\xi\left( {\xi - 1} \right)}^{2}\left( {s_{j + 1} - s_{j}} \right)}}{{f_{3j}(s)} = {\left( {3 - {2\xi}} \right)\xi^{2}}}{{f_{4j}(s)} = {{\xi^{2}\left( {\xi - 1} \right)}\left( {s_{j + 1} - s_{j}} \right)}}{\xi = \frac{s - s_{j}}{s_{j + 1} - s_{j}}}} & (32)\end{matrix}$The cubic spline functions defined in equation (32) are not the onlypossible choices. An alternate formulation that has direct connection todrillstring solutions is the tension spline:

$\quad\begin{matrix}{{{f_{1j}(\xi)} = {1 + \frac{\left. \left\lbrack {{\cosh(\lambda)} - 1} \right) \right\rbrack\left\lbrack {1 - {\cosh({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sinh(\lambda)}} + {2\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}} - \frac{{\sinh(\lambda)}\left\lbrack {{\lambda\xi} - {\sinh({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sinh(\lambda)}} + {2\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}}}{{f_{2j}(\xi)} = {\left\{ {\xi - \frac{\left. \left\lbrack {{\sinh(\lambda)} - {\lambda\;{\cosh(\lambda)}}} \right) \right\rbrack\left\lbrack {1 - {\cosh({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sinh(\lambda)}} + {2{\lambda\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}} - \frac{\left\lbrack {{\lambda\;\sinh\;(\lambda)} + 1 - {\cosh(\lambda)}} \right\rbrack\left\lbrack {{\lambda\xi} - {\sinh({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sinh(\lambda)}} + {2{\lambda\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}}} \right\}\left( {s_{j + 1} - s_{j}} \right)}}{{f_{3j}(\xi)} = {\frac{\left. \left\lbrack {1 - {\cosh(\lambda)}} \right) \right\rbrack\left\lbrack {1 - {\cosh({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sinh(\lambda)}} + {2\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}} + \frac{{\sinh(\lambda)}\left\lbrack {{\lambda\xi} - {\sinh({\lambda\xi})}} \right\rbrack}{{{\lambda sinh}(\lambda)} + {2\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}}}{{f_{4j}(\xi)} = {\left\{ {\frac{\left. \left\lbrack {{\sinh(\lambda)} - \lambda} \right) \right\rbrack\left\lbrack {1 - {\cosh({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sinh(\lambda)}} + {2{\lambda\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}} + \frac{\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack\left\lbrack {{\lambda\xi} - {\sinh({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sinh(\lambda)}} + {2{\lambda\left\lbrack {1 - {\cosh(\lambda)}} \right\rbrack}}}} \right\}\left( {s_{j + 1} - s_{j}} \right)}}\mspace{79mu}{\xi = \frac{s - s_{j}}{s_{j + 1} - s_{j}}}} & (33)\end{matrix}$where λ is a parameter to be determined. For beam-column solutions,

$\begin{matrix}{{\lambda = {\Delta\; s\sqrt{\frac{\overset{\sim}{F}}{EI}}}}{\overset{\sim}{F} = {{F - {{EI}\;\kappa^{2}}} > 0}}} & (34)\end{matrix}$A similar solution for strings in compression is:

$\begin{matrix}{{{f_{1j}(\xi)} = {1 - \frac{\left. \left\lbrack {{\cos(\lambda)} - 1} \right) \right\rbrack\left\lbrack {1 - {\cos({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sin(\lambda)}} - {2\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}} - \frac{{\sin(\lambda)}\left\lbrack {{\lambda\xi} - {\sin({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sin(\lambda)}} - {2\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}}}{{f_{2j}(\xi)} = {\left\{ {\xi + \frac{\left. \left\lbrack {{\sin(\lambda)} - {\lambda\;{\cos(\lambda)}}} \right) \right\rbrack\left\lbrack {1 - {\cos({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sin(\lambda)}} - {2{\lambda\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}} - \frac{\left\lbrack {{\lambda\;\sin\;(\lambda)} - 1 + {\cos(\lambda)}} \right\rbrack\left\lbrack {{\lambda\xi} - {\sin({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sin(\lambda)}} - {2{\lambda\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}}} \right\}\left( {s_{j + 1} - s_{j}} \right)}}{{f_{3j}(\xi)} = {\frac{\left. \left\lbrack {{\cos(\lambda)} - 1} \right) \right\rbrack\left\lbrack {1 - {\cos({\lambda\xi})}} \right\rbrack}{{\lambda\;{\sin(\lambda)}} - {2\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}} + \frac{{\sin(\lambda)}\left\lbrack {{\lambda\xi} - {\sin({\lambda\xi})}} \right\rbrack}{{{\lambda sin}(\lambda)} - {2\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}}}{{f_{4j}(\xi)} = {\left\{ {{- \frac{\left. \left\lbrack {{\sin(\lambda)} - \lambda} \right) \right\rbrack\left\lbrack {1 - {\cos({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sin(\lambda)}} - {2{\lambda\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}}} + \frac{{\lambda\left\lbrack {{\cos(\lambda)} - 1} \right\rbrack}\left\lbrack {{\lambda\xi} - {\sin({\lambda\xi})}} \right\rbrack}{{\lambda^{2}{\sin(\lambda)}} - {2{\lambda\left\lbrack {1 - {\cos(\lambda)}} \right\rbrack}}}} \right\}\left( {s_{j + 1} - s_{j}} \right)}}\mspace{79mu}{\xi = \frac{s - s_{j}}{s_{j + 1} - s_{j}}}} & (35)\end{matrix}$where λ is a parameter to be determined. For beam-column solutions,

$\begin{matrix}{{\lambda = {\Delta\; s\sqrt{\frac{- \overset{\sim}{F}}{EI}}}}{\overset{\sim}{F} = {{F - {{EI}\;\kappa^{2}}} < 0}}} & (36)\end{matrix}$

One problem is that the λ coefficients are functions of the axial force,which are not known until the torque-drag equations are solved. Inpractice, λ tends to be small, so that the solution approximates a cubicequation. The cubic interpolation can be used to approximate thetrajectory, and to solve the torque-drag problem. The torque-dragsolution can then be used to refine the trajectory, iterating ifnecessary.

A simple comparison of the wellbore trajectory model of the presentinvention, also referred to as a spline model, and the standard minimumcurvature model with three analytic wellbore trajectories (circular-arc,catenary, helix) is illustrated in FIGS. 2-4, respectively. Thecomparisons of the displacements illustrated in FIGS. 2-4 demonstratethat the minimum curvature model and the spline model match the analyticwellbore trajectory in FIG. 2 (circular-arc), the analytic wellboretrajectory in FIG. 3 (catenary) and the analytic wellbore trajectory inFIG. 4 (helix). Only one displacement is shown for the helix, but isrepresentative of the other displacements. The spline model was alsoused to calculate the rate of change of curvature for the analyticwellbore trajectory illustrated in FIG. 5 (catenary), and the geometrictorsion for the analytic wellbore trajectory illustrated in FIG. 6(helix). Despite the results of the simple comparison illustrated inFIGS. 2-4, the results illustrated by the comparisons in FIGS. 5-6demonstrate the deficiencies of the minimum curvature model whencalculating the curvature rate of change for the catenary wellboretrajectory illustrated in FIG. 5 or when calculating the geometrictorsion for the helix wellbore trajectory illustrated in FIG. 6. Theminimum curvature model predicts zero for both quantities compared inFIGS. 5-6, which cannot be plotted. The spline model, however,determines both quantities accurately, although there is some end effectapparent in the geometric torsion calculation. Additional advantagesattributed to the present invention (spline model) are demonstrated bythe following examples.

Torque-Drag Calculations

Torque-drag calculations were made using a comprehensive torque-dragmodel well known in the art. Similarly, the equilibrium equations wereintegrated using a method well known in the art. Otherwise, the onlydifference in the solutions is the choice of the trajectory model.

Example 1

In this example, the drag and torque properties of an idealized wellplan are based on Well 3 described in Society of Petroleum Engineersarticle “Designing Well Paths to Reduce Drag and Torque” by Sheppard, M.C., Wick, C. and Burgess, T. M. Referring now to FIG. 7, the fixedpoints on the model trajectory are as follows: i) the well is consideredto be drilled vertically to a KOP at a depth of 2,400 ft.; ii) theinclination angle then builds at a rate of 5°/100 ft; and iii) thetarget location is considered to be at a vertical depth of 9,000 ft anddisplaced horizontally from the rig location by 6,000 ft. Drilled as aconventional build-tangent well, this would correspond to a 44.5° welldeviation. The model drillstring was configured with 372 feet of 6½ inchdrill collar (99.55 lbf/ft.) and 840 ft of 5 inch heavyweight pipe(50.53 lbf./ft.) with 5 inch drillpipe (20.5 lbf./ft.) to the surface. Amud weight of 9.8 lbm/gal was used. In this example, a value of 0.4 waschosen for the coefficient of friction to simulate severe conditions.Torque-loss calculations were made with an assumed WOB of 38,000 lbf.and with an assumed surface torque of 24,500 ft.-lbf.

Hook load calculated for zero friction was 192202 lbf. for thecircular-arc calculation, and 192164 lbf. for the spline model, whichcompare to a spreadsheet calculation of 192203 lbf. The slightdifference (38 lbf.) is due to the spline taking on a slightly differentshape (due to smoothness requirements) from thestraight-line/circular-arc shapes specified, which the minimum curvaturemodel exactly duplicated. Other than the slight difference in the splinetrajectory, all other aspects of the axial force calculations areidentical between the two models. Tripping out, with a frictioncoefficient of 0.4, the hook load was 313474 lbf. for the circular-arcmodel and 319633 lbf. for the spline model, for a difference of 6159lbf. If calculations are from the zero friction base line, thisrepresents a difference of 5% in the axial force loading. With a surfacetorque of 24,500 ft-lbs., the torque at the bit was 3333 ft-lbs. for theminimum curvature model and 2528 ft-lbs. for the spline model. Thisrepresents a 4% difference in the distributed torque between the twomodels. The bending moments for the drillstring through the buildsection are illustrated in FIG. 8. Notably, the minimum curvature doesgive a lower bending moment than the spline, but that the spline resultsare much smoother.

Since this case has a relatively mild build rate, and since the buildsection was only about 8% of the total well depth, it would be expectedthat a relatively small effect from the spline formulation would beseen. Because the classic torque-drag analysis has historically givengood results, the agreement of the two models for this case verifiesthat the overall formulation is correct.

Example 2

For a more demanding example, the short-radius wellbore described in theSociety of Petroleum Engineers article “Short Radius TTRD Well with RigAssisted Snubbing on the Veslefrikk Field” by Grinde, Jan, and Haugland,Torstein was used. Referring now to FIGS. 9A and 9B, the vertical andhorizontal views of the end of the wellpath are illustrated,respectively. The build rate for this example was 42°/30 m, roughly tentimes the build rate of the first case in Example 1. As illustrated inFIG. 10, some of the contact force is neglected by neglecting thebending moment since the contact force for the spline model at the endof the build is four times that of the minimum curvature model. In FIG.11, the bending moment for this example is illustrated. The minimumcurvature model still provides a lower bending moment than the splinemodel, but the spline results are still much smoother.

Referring now to FIG. 12, flow diagram illustrates one embodiment of amethod 1200 for implementing the present invention.

In step 1202, the survey data is obtained for each survey point (_(j)).

In step 1204, a tangent vector ({right arrow over (t)}_(j)) iscalculated at each survey point using the survey data at each respectivesurvey point.

In step 1206, a special normal vector (ñ_(j)) and a special binormalvector ({tilde over (b)}_(j)) are calculated at each survey point.

In step 1208, a block tridiagonal matrix is calculated using the tangentvector, the special normal vector and the special binormal vector ateach respective survey point.

In step 1210, a coefficient (α_(j)) is calculated at each survey pointin the direction of the special normal vector at the respective surveypoint and another coefficient (β_(j)) is calculated at each survey pointin the direction of the special binormal vector at the respective surveypoint using the block tridiagonal matrix.

In step 1212, a wellbore curvature (κ_(j)) and a normal vector ({rightarrow over (n)}_(j)) are calculated at each survey point using a firstderivative of the tangent vector, the coefficient and the anothercoefficient at each respective survey point.

In step 1214, a tangent vector interpolation function ({right arrow over(t)}_(j)(s)) is calculated for each interval between survey points usingthe wellbore curvature, the tangent vector and the normal vector at eachrespective survey point.

In step 1216, the wellbore trajectory is determined using each tangentvector interpolation function in a torque-drag drillstring model.

While the present invention has been described in connection withpresently preferred embodiments, it will be understood by those skilledin the art that it is not intended to limit the invention to thoseembodiments. The present invention, for example, may also be applied tomodel other tubular trajectories, which are common in chemical plantsand manufacturing facilities. It is therefore, contemplated that variousalternative embodiments and modifications may be made to the disclosedembodiments without departing from the spirit and scope of the inventiondefined by the appended claims and equivalents thereof.

The invention claimed is:
 1. A computer implemented method for modelinga wellbore trajectory, comprising: calculating a tangent vectorinterpolation function for each interval between two or more surveypoints within a wellbore; and determining the wellbore trajectory usinga computer processor and each tangent vector interpolation function in atorque-drag drillstring model.
 2. The method of claim 1, furthercomprising: calculating a tangent vector at each survey point usingsurvey data at each respective survey point, the tangent vector beingused to calculate the tangent vector interpolation function.
 3. Themethod of claim 2, wherein the survey data comprises an inclinationangle, an azimuth angle and a measured depth at each survey point. 4.The method of claim 1, further comprising: calculating a wellborecurvature at each survey point using a first derivative of the tangentvector, a coefficient and another coefficient at each respective surveypoint, the wellbore curvature being used to calculate the tangent vectorinterpolation function; and calculating a normal vector at each surveypoint using the first derivative of the tangent vector, the coefficientand the another coefficient at each respective survey point, the normalvector being used to calculate the tangent vector interpolationfunction.
 5. The method of claim 4, wherein the first derivative of thetangent vector is continuous at each survey point.
 6. The method ofclaim 4, further comprising: calculating the coefficient at each surveypoint in a direction of a special normal vector at the respective surveypoint using a block tridiagonal matrix; and calculating the anothercoefficient at each survey point in a direction of a special binormalvector at the respective survey point using the block tridiagonalmatrix.
 7. The method of claim 6, further comprising: calculating theblock tridiagonal matrix using the tangent vector, the special normalvector, and the special binormal vector at each respective survey point.8. The method of claim 6, further comprising: calculating the specialnormal vector at each survey point; and calculating the special binormalvector at each survey point.
 9. The method of claim 1, furthercomprising: calculating a torque-drag drillstring solution using thewellbore trajectory.
 10. The method of claim 9, further comprising:refining the wellbore trajectory using the torque-drag drillstringsolution.
 11. A non-transitory computer readable medium having computerexecutable instructions for modeling a wellbore trajectory, theinstructions being executable to implement: calculating a tangent vectorinterpolation function for each interval between two or more surveypoints within a wellbore; and determining the wellbore trajectory usingeach tangent vector interpolation function in a torque-drag drillstringmodel.
 12. The computer readable medium of claim 11, further comprising:calculating a tangent vector at each survey point using survey data ateach respective survey point, the tangent vector being used to calculatethe tangent vector interpolation function.
 13. The computer readablemedium of claim 12, wherein the survey data comprises an inclinationangle, an azimuth angle and a measured depth at each survey point. 14.The computer readable medium of claim 11, further comprising:calculating a wellbore curvature at each survey point using a firstderivative of the tangent vector, a coefficient and another coefficientat each respective survey point, the wellbore curvature being used tocalculate the tangent vector interpolation function; and calculating anormal vector at each survey point using the first derivative of thetangent vector, the coefficient and the another coefficient at eachrespective survey point, the normal vector being used to calculate thetangent vector interpolation function.
 15. The computer readable mediumof claim 14, wherein the first derivative of the tangent vector iscontinuous at each survey point.
 16. The computer readable medium ofclaim 14, further comprising: calculating the coefficient at each surveypoint in a direction of a special normal vector at the respective surveypoint using a block tridiagonal matrix; and calculating the anothercoefficient at each survey point in a direction of a special binormalvector at the respective survey point using the block tridiagonalmatrix.
 17. The computer readable medium of claim 16, furthercomprising: calculating the block tridiagonal matrix using the tangentvector, the special normal vector, and the special binormal vector ateach respective survey point.
 18. The computer readable medium of claim16, further comprising: calculating the special normal vector at eachsurvey point; and calculating the special binormal vector at each surveypoint.
 19. The computer readable medium of claim 11, further comprising:calculating a torque-drag drillstring solution using the wellboretrajectory.
 20. The computer readable medium of claim 19, furthercomprising: refining the wellbore trajectory using the torque-dragdrillstring solution.